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Using a combination of Monte Carlo techniques, we locate the liquid-vapor critical point of 
adhesive hard spheres. We find that the critical point lies deep inside the gel region of the phase 
diagram. The (reduced) critical temperature and density are Tc — 0.1133 ± 0.0005 and Pc = 0.508 ± 
0.01. We compare these results with the available theoretical predictions. Using a finite-size scaling 
analysis, we verify that the critical behavior of the adhesive hard sphere model is consistent with 
that of the 3D Ising universality class. 
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The structure of a simple liquid is well described by 
that of a system of hard spheres at the same effective 
density. To a good approximation, the effect of attractive 
interactions on the liquid structure can be ignored. This 
feature of simple liquids is implicit in the Van der Waals 
theory for the liquid-vapor transition, and has been made 
explicit in the highly successful thermodynamic pertur- 
bation theories for simple liquids 1]. The perturbation 
approach becomes exact as the range of the attractive 
interaction tends to infinity while its integrated strength 
remains constant We refer to this limit as the 'Van 
der Waals' (VDW) Hmit 0. Conversely, as the attrac- 
tive forces become shorter-ranged and stronger, the per- 
turbation approach is likely to break down. Fluids with 
strong, short-ranged attraction (so-called 'energetic' flu- 
ids 4!|) are of growing importance in the area of complex 
liquids. For example, short-range attractions are thought 
to be responsible for the transition from a 'repulsive' to 
an 'attractive' glass which has recently been observed 
experimentally in PMMA (polymethylmethacrylate) dis- 
persions 0- 

In this Letter, we consider a model system that can be 
considered as the prototypical energetic fluid: a fluid of 
adhesive hard spheres (AHS). Introduced in 1968 0, the 
AHS model is a reference system for particles with short 
range attractions. The pair potential consists of an im- 
penetrable core plus a surface adhesion term that favors 
configurations where spheres are in contact. At larger 
separations, there is no interaction. The AHS model can 
be considered as the 'anti Van der Waals' limit. 

Baxter showed '7] that the Percus-Yevick (FY) equa- 
tion can be solved analytically for adhesive hard spheres. 
In fact, Baxter's solution is often used to to analyze ex- 
perimental results for systems as diverse as silica suspen- 
sions Q, copolymer micelles lit, and the fluid phase of 
lysozyme jlflj . 

One important feature of the AHS model is that its 
phase diagram contains a liquid-vapor coexistence re- 
gion 01 . The FY equation offers different routes to 
estimate the location of the liquid-vapor critical point. 
However, the 'compressibility route' l7l| and the 'energy 



route' [13 lead to estimates for the critical temperature 
that differ by some 20%, while the estimates for the crit- 
ical density differ by almost a factor of three. For the 
analysis of experimental data, it is important to know 
the location of the critical point more accurately. The 
reason is that, upon cooling, adhesive hard spheres can 
percolate to form a gel. Information about the location of 
this percolation curve is available from simulation 13, l"3| 
and from analytical estimates [llij . It is not clear, though, 
how percolation interferes with the fluid-fluid phase sep- 
aration. If Baxter's estimate of the critical point were 
right, the critical point would be on, or near, t he p erco- 
lation curve. However, if the estimate of Ref. 12] were 
closer to reality, then the fluid-fluid critical point would 
be deep inside the gel phase. Clearly, this difference 
has implications for the possibility to observe the fluid- 
fluid critical point in systems with short-ranged attrac- 
tion. Also for the (attractive) glass transition, it makes a 
considerable difference whether or not the transition line 
runs close to a liquid- vapor critical point. 

In spite of its obvious importance as a reference sys- 
tem for 'energetic' (as contrasted to 'entropic') complex 
liquids Q, there exist, to our knowledge, no accurate nu- 
merical estimates of the critical point of the AHS model. 
This is not surprising, as computer simulations of the low- 
temperature AHS model are notoriously difficult, due its 
propensity to form large, even percolating, clusters. In 
this Letter, we report a Monte Carlo (MC) simulation 
study that allows us to locate the AHS critical point. 

The AHS interaction is derived from a square well po- 
tential by taking a limit in which the well becomes in- 
finitesimally narrow at the same time as becoming in- 
finitely deep in such a way that the integrated Boltzmann 
weight of bound configurations remains finite ■ The in- 
teraction is most easily defined by the expression for the 
Boltzmann factor as a function of pair separation, r: 

exp[-C/(r)/fcr] = e(r - a) + -^S{r - a). (1) 

12r 

In Eq. U{r) and kT are the formal pair potential and 
thermal energy, respectively, while a is the hard core di- 
ameter, which we henceforth use as the unit of length. 
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The step function 8 accounts for the hard sphere repul- 
sion, and the Dirac delta function introduces the surface 
adhesion. The parameter r determines the strength of 
the adhesion and can be interpreted as an effective tem- 
perature, or an 'inverse stickiness parameter.' 

To compare the properties of the AHS and VDW mod- 
els, we need a measure for the strength of the attractive 
interactions that is meaningful in both limits. This is 
best achieved by comparing the reduced second virial 
coefficients: B2 = 82/8^'^, where B^^ is the second 
virial coefficients of hard spheres. For the AHS model, 
B2 = I — 1/(4t), and we use this expression to define the 
T-parameter for the VDW model. In the VDW limit, the 
free-energy density is given by /vdw(p) = /hs(p) - ap^ 
where /hs is the free energy density of a system of hard 
spheres at density p, while the constant a measures the 
strength of the attractive forces. Hence, in the VDW 
limit, T = kTB^^ /(Aa) is proportional to the tempera- 
ture. We henceforth refer to r as the temperature pa- 
rameter for both the VDW and AHS systems. 

In the AHS model, particles may stick and form clus- 
ters. At sufficiently low r, clusters may percolate (span 
the system), mimicking the infinite clusters that form 
during gelation in the real system. Percolating clusters 
pose serious problems for any simulation scheme that 
employs volume changing moves (e.g. constant pressure 
MC), since all such moves are rejected as soon as perco- 
lating clusters appear. 

To study the phase behavior of adhesive hard spheres, 
we therefore used Grand Canonical MC (GCMC) sim- 
ulation. This technique can be used beyond the perco- 
lation threshold. Yet, even with GCMC, equilibration 
of the system at low temperatures is prohibitively slow. 
We therefore combine GCMC with parallel tempering 
to speed up equilibration, and with multiple histogram 
reweighting, to make optimal use of the available simula- 
tion data. All simulation techniques had to be specifically 
adapted for the AHS model. 

The fact that the attractive term in Eq. is infinites- 
imally narrow means that, in a simulation, the chance of 
generating a bound configuration by random displace- 
ment of a particle is vanishingiy small. Conversely, the 
probability of breaking such a bond, once formed, is also 
negligible. Consequently, conventional Metropolis sam- 
pling cannot be applied to the AHS model. 

We therefore employ a modification of the AHS-MC 
schemes described in Ref. jlSj, where single particle dis- 
placements explicitly make and break up to three con- 
tacts with the test particle simultaneously. States with 
higher coordination numbers can be reached indirectly 
through suitable combinations of moves. 

In our GCMC simulations, we employ particle inser- 
tion and removal steps with equal probability that to- 
gether constitute 45% of MC steps. We only insert and 
remove particles that are not bound to others. 

To speed up equilibration, we perform cluster transla- 



tion moves with probability 5%. A particle is chosen at 
random, and the cluster to which it belongs is translated 
by a random amount in each Cartesian direction up to 
a maximum that is inversely proportional to the number 
of particles in the cluster. 

To overcome the slow equilibration of large clusters at 
low T or high density, we have used the parallel tem- 
pering scheme of Geyer Parallel tempering can be 
performed with replicas at different temperatures, differ- 
ent chemical potentials, or combinations of both. We 
have found it most convenient to gather statistics across 
the full density range for one temperature at a time, and 
therefore chose to run a hierarchy of chemical potentials, 
/i, with the same value of r. Each replica attempts a 
configuration exchange with one of its neighbors in fi (al- 
ternating between the higher and lower neighbor) every 
200 MC steps. The acceptance probability for an ex- 
change between replicas i and j is min[l, {zi/ Zj)^^'^^], 
where Ni is the number of particles in the current con- 
figuration of replica i, and the activity z is related to the 
chemical potential by z = A~'^ exp(/i/A:T), with A the 
thermal de Broglie wavelength. 

The final computational tool needed to analyze the 
data is multiple histogram reweighting 0. The proba- 
bility of observing the system in a state with N particles 
and B binary contacts between the particles at tempera- 
ture T and activity z can be exactly decomposed into the 
form 

PNB{T,z) = nNBT~''z''/E{T,z), (2) 

where ^Inb is an effective density of states, and S(r, z) 
is the grand canonical partition function. Knowledge of 
f^Afs up to a multiplicative constant is sufficient to cal- 
culate the distribution oi N or B at any set of conditions 
(t, z). flNB can be obtained over a wide range of {N, B) 
by combining overlapping two-dimensional histograms of 
N and B collected from simulations at different r and z. 

Simulations at sufficiently low r show a range of activ- 
ities in which the density distribution is bimodal, indi- 
cating the coexistence of a high- and a low-density fluid 
phase. We studied this coexistence at four different sys- 
tem sizes, labeled by the length L of the cubic simulation 
box. For L = 5, 6 and 8, the simulations consisted of 10^ 
equilibration blocks and 10^ acquisition blocks, where a 
block is ja^ MC trial moves, each of which may be a 
particle displacement, particle insertion/removal, cluster 
translation, or replica exchange. At each temperature, 
seven or eight parallel replicas were used, with the range 
of activities chosen to span the coexisting densities. At 
L — 10, the numbers of equilibration and acquisition 
blocks had to be lowered to 2.5 x 10^ and 2.5 x 10^, 
respectively, to obtain results in a reasonable computa- 
tional time. 

The (iV, i?) histograms for different t and z were com- 
bined for each system size separately. The coexisting den- 
sities were then obtained as a function of temperature by 
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TABLE I: Size-dependent properties for critical point deter- 
mination. Ta is the apparent critical temperature, at which 
the distribution of A4 (with field-mixing parameter s) col- 
lapses onto the universal form p*{x). The remaining three 
columns refer to properties at the proposed critical tempera- 
ture of t"™ — 0.1133: pc is the mean density, Zc is the activity, 
and aL scales p(ai,A^) to have unit variance. 
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FIG. 1: Coexistence curve of the AHS model from GCMC 
using four box lengths, L. Solid symbols denote points ac- 
tually simulated; open symbols were obtained by histogram 
reweighting. Error bars are omitted for clarity. The cross 
denotes the estimated critical point at r = 0.1133, p — 0.508. 



reweighting the histograms to find the activity at which 
the two density peaks have equal height. Figure ^ shows 
the resuhing coexistence curves. Strong finite size effects 
are immediately visible. The curves obtained from differ- 
ent system sizes coincide at sufficiently low temperature, 
but deviate significantly from each other as the critical 
region is approached. 

The size dependence of the coexistence distributions 
can be used to locate the critical point more accurately 
using the approach developed by Bruce and Wilding 
[lilfl^ . This method uses the fact that, within a univer- 
sality class, the critical distribution of the order parame- 
ter is invariant up to a rescaling of the order parameter. 
As the AHS interactions are short-ranged, the fluid-fluid 
critical point of this model is expected to belong to the 
3D Ising universality class. The critical distribution is 
known to high precision from studies of lattice systems, 
and accurate analytic fits are available [20| . 

Due to the absence of particle-hole symmetry in off- 
lattice models, the distribution of the density in the AHS 
model is not quite symmetric about its mean. Symmetry 
can be restored by accounting for the mixed character 
of the scaling fields. The appropriate order parameter is 
then no longer the pure particle density p — N/L^, but 
includes a contribution from the energy density, which 
in the AHS model can be taken as u = —B/L^. The 
particle and energy density operators are replaced by the 
linear combinations |l8j 
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where s and r are system-dependent field mixing param- 
eters that are identically zero for models with Ising sym- 
metry. Precisely at criticality, the distribution of A4 in 
a system of sufficiently large linear length L takes on 



the universal form pl{M) — P*{aL[J^ — -^c]), where 
Aic — i-M) evaluated at the critical temperature, and 
aL oc L^/" with (3 and v critical exponents. 

Tabled lists the apparent critical temperature Ta and 
mixing parameter s at which the distribution of M most 
closely matches p*{x) for the four system sizes studied. 
The fact that s is small indicates that the distribution of p 
itself is nearly symmetric. The apparent size-dependence 
of s is probably — at least partly — due to the fact that, 
despite the long simulations employed, statistical fluctu- 
ations in the histograms are still signiflcant. For small 
L, Ta is expected to show some size dependence due to 
finite-size corrections to scaling |0. However, the differ- 
ences in Ta shown in Table d^^re not monotonic, and are 
again likely to be due to limited statistics. Rather than 
attempting to extrapolate to the infinite system limit, we 
therefore quote the critical parameters that most closely 
reproduce the universal critical form, with conservative 
error bars: 



0.1133 ±0.0005 



0.508 ±0.01. 



Figure |21 shows the distribution of the order parameter 
at the proposed r^'™ (not the individual apparent critical 
temperatures Ta) for all four system sizes. The collapse 
of the data onto the Ising distribution shows that the 
transition is consistent with this universality class. 

The size dependence of the scaling factor in Table HI 
permits an estimate of the ratio of the critical exponents 
(3 and v. A straightforward fit to the asymptotic form 
ql (X Ll^/" yields P/v = 0.50 ± 0.04, which is clearly 
compatible with the Ising value (i/v — 0.52. 

We can now compare our numerical results with var- 
ious theoretical predictions. Below, we list the esti- 
mates based on the PY compressibility and energy 
[l^ routes, and the estimate based on the VDW (mean- 
field) expression [2ll |: 
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FIG. 2: Distribution of the order parameter at the proposed 
critical temperature from GCMC at four different box lengths, 
L. Each curve has been scaled to have unit variance. The 
solid curve is the 3D Ising critical distribution l23|. 



lation lines clearly show that the critical point of the 
fluid-fluid transition lies well within the percolated part 
of the phase diagram. 

The phase diagram reported in Fig. |31 is likely to be 
representative of that of real colloidal systems with short- 
ranged attraction at the same value of r. If so, our find- 
ings suggest that the attractive glass transition 0, |^ is 
relatively far removed from the fluid-fluid critical point. 
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FIG. 3: Solid lines: coexistence curves for the VDW limit, the 
AHS PY compressibility and energy routes, and AHS simula- 
tion with L = 8. For the latter, the line merely guides the eye 
between the data points. Dashed lines: percolation threshold 
from PY theory and simulation (line with circles). 



The coexistence curves for the two PY routes can be ob- 
tained numerically using the analytical expressions 
for the chemical potential and pressure. These theoret- 
ical curves, together with the mean-field result and the 
simulation data at i = 8 are shown in Fig. |31 The sim- 
ulation results lie between the two PY predictions, but 
are clearly closer to those of the energy equation. 

Also plotted in Fig.|31is the percolation threshold. The 
theoretical result corresponds to the PY estimate of 
the points where the cluster size diverges. In simulations, 
we define the percolation threshold as the locus of points 
were the probability of observing a percolating cluster 
in a canonical simulation is 50%. This definition is rel- 
atively size independent, and our results with N — 500 
particles are consistent with the more elaborate analysis 
of Lee 01 . Both the theoretical and simulated perco- 
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